11.04.2022

Linear modeling

model predictors
Regression continuous
ANOVA categorical
ANCOVA categorical & continuous
Mixed-effects models categorical (& continuous) with random effects


Examples with lm() models.
But glm() (generalized linear models: Poisson, logistic, binomial, …) works similarly

Linear modeling

summary table after fitting a model with lm()

summary.lm() or summary() summary.aov() or anova()
Regression t-tests for effect sizes F-tests for stepwise nested models
ANOVA t-tests for dummy-coded effects (levels) F-tests for stepwise nested models
ANCOVA t-tests for effect sizes and dummy-coded effects (levels) F-tests for stepwise nested models

Overview

  1. Hypothesis testing
  2. Regression
  3. ANOVA
  4. ANCOVA
  5. Mixed models

Overview

  1. Hypothesis testing
  2. Regression
  3. ANOVA
  4. ANCOVA
  5. Mixed models

Hypothesis testing

  • assumptions: random variable \(X\sim\text{normal}(\mu,\sigma)\)
  • data \(x_1,\dots,x_n\)
  • H0: Null hypothesis with random variable \(X\) - \(\mu=\text{E}(X)=0\)
  • H1: Alternative hypothesis - e.g. \(\mu>0\)
  • Test statistic \(T(X)\) transformation
    distribution of \(T\) under H0 is known - \(T(X)=\frac{\bar X-\mu}{\sigma}=\frac{1}{\sigma}\frac{1}{n}\sum_{i=1}^n{x_i}\sim\text{normal}(0,1)\)

Hypothesis testing

  • define type I error \(\alpha\) (e.g. \(\alpha=0.05\)):
    chance of rejecting H0 (\(\mu=0\)) although it is true
  • check if \(T(\text{data})\) is outside \((1-\alpha)\) quantile
    âž” reject H0!
  • P-value: \(P(T>T(\text{data}))\)

Z-test

\(X_i\) have normal distribution

\(Z=\sum_i X_i\) sum has normal distribution

e.g. Z-test for mean (standard deviation \(\sigma\) is known)

Chi-squared-test

\(X_i\) have standard normal distribution

\(Z=\sum_iX_i^2\) sum of squares has \(\chi^2_n\) distribution

e.g. chisq-test for independence

F-test

\(X\sim\chi^2_n\), \(Y\sim\chi^2_m\)

\(Z=\displaystyle\frac{X/n}{Y/m}\) ratio has \(F_{n,m}\) distribution

e.g. F-test for regression, ANOVA, for comparing variances

t-test

\(X\) normal, \(Y\sim\chi_n^2\)

\(Z=\displaystyle\frac{X}{\sqrt{Y/n}}\) has \(t_n\) distribution

e.g. t-test for regression coefficients, 2-sample t-test for means

Palmer penguins

Artwork by @allison_horst

Palmer penguins

library(palmerpenguins)
data=as.data.frame(penguins)
str(data)
## 'data.frame':    344 obs. of  8 variables:
##  $ species          : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ island           : Factor w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ bill_length_mm   : num  39.1 39.5 40.3 NA 36.7 39.3 38.9 39.2 34.1 42 ...
##  $ bill_depth_mm    : num  18.7 17.4 18 NA 19.3 20.6 17.8 19.6 18.1 20.2 ...
##  $ flipper_length_mm: int  181 186 195 NA 193 190 181 195 193 190 ...
##  $ body_mass_g      : int  3750 3800 3250 NA 3450 3650 3625 4675 3475 4250 ...
##  $ sex              : Factor w/ 2 levels "female","male": 2 1 1 NA 1 2 1 2 NA NA ...
##  $ year             : int  2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...

Example

2 sample t-test

Data in 2 groups: \(x_1\dots x_n\)
                               \(y_1\dots y_m\)

Assumption: \(X\sim\text{normal}(\mu_X,\sigma_X)\)
                        \(Y\sim\text{normal}(\mu_Y,\sigma_Y)\)

H0: \(\mu_X=\mu_Y\) (means are identical)

H1: \(\mu_X\neq\mu_Y\)

under H0: \(T=( \bar X - \bar Y ) / \sqrt{\frac{s^2_X}{n}+\frac{s^2_Y}{m}}\) has t-distribution

Example

x = subset(data, species=="Chinstrap")$flipper_length_mm
y = subset(data, species=="Adelie")$flipper_length_mm
boxplot(x,y)

Example

t.test(x,y)
## 
##  Welch Two Sample t-test
## 
## data:  x and y
## t = 5.7804, df = 119.68, p-value = 6.049e-08
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  3.859244 7.880530
## sample estimates:
## mean of x mean of y 
##  195.8235  189.9536

âž” reject H0 (\(P<0.05\)) that the means are identical

Example

## t = 5.7804, df = 119.68, p-value = 6.049e-08

curve(dt(x, df=119.68), from=0, to=6, lwd=3, col=4)
abline(v=qt(p=0.95, df=119.68), lty=2  )
points(5.7804, 0, col=2, pch=16, cex=1.5)

Overview

  1. Hypothesis testing
  2. Regression
  3. ANOVA
  4. ANCOVA
  5. Mixed models

Overview

summary table after fitting a model with lm()

summary.lm() or summary() summary.aov() or anova()
Regression t-tests for effect sizes F-tests for stepwise nested models
ANOVA t-tests for dummy-coded effects (levels) F-tests for stepwise nested models
ANCOVA t-tests for effect sizes and dummy-coded effects (levels) F-tests for stepwise nested models

Sums of squares

\(SST=\sum_i({y}_i-\bar{y})^2\)     total
\(SSE=\sum_i(\widehat{y}_i-\bar{y})^2\)    explained
\(SSR=\sum_i({y}_i-\widehat{y}_i)^2\)   residual

\(SST=SSE+SSR\)

Goodness of fit

variation explained by the model

\(R^2=\frac{SSE}{SST}=1-\frac{SSR}{SST}\)

F-test comparing two models

  • two nested models, e.g.    m1: y~x1*x2
                                                   m2: y~x1+x2
  • m1 better than m2?
    residuals: \(SSR_1<SSR_2\)?
  • H0: \(SSR_1=SSR_2\)
  • test statistic \(\frac{(SSR_2-SSR_1)/(k_1-k_2)}{SSR_1/(n-k_2-1)}\) has \(F_{k_1-k_2,n-k_2-1}\) distr.
    (ratio of two \(\chi^2\) distributions)
  • If \(P<0.05\) âž” reject H0
  • H1: m1 has better fit

F-test comparing two models

use anova( , )

data = subset(data, species=="Chinstrap")
mod1 = lm(flipper_length_mm ~ body_mass_g * bill_length_mm, data)
mod2 = lm(flipper_length_mm ~ body_mass_g + bill_length_mm, data)
anova(mod1, mod2)
## Analysis of Variance Table
## 
## Model 1: flipper_length_mm ~ body_mass_g * bill_length_mm
## Model 2: flipper_length_mm ~ body_mass_g + bill_length_mm
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     64 1901.2                           
## 2     65 1911.8 -1   -10.585 0.3563 0.5527
  • H0: mod1 and mod2 perform equally
  • \(P>0.05\), can’t reject H0
  • choose simpler model mod2

F-test comparing two models

## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 64 1901.2
## 2 65 1911.8 -1 -10.585 0.3563 0.5527

curve(df(x, df1=1, df2=64), from=0, to=4, col=4, lwd=2)
abline(v=qf(p=0.95, df1=1, df2=64), lty=2)
points(0.3563, 0, col=2, pch=16, cex=1.5)

Overall F-test for regression

  • test your regression model     y~1+x1+...+xk
    against intercept only model  y~1
  • H0: regression model is only as good as intercept model
           same as \(SSR=SST\)
           same as \(R^2=0\)
           same as \(\beta_1=\ldots=\beta_k=0\)
  • test statistic \(\frac{SSE/k}{SSR/(n-k-1)}\) has \(F_{k,n-k-1}\) distribution
  • If \(P<0.05\) âž” reject H0
  • H1: at least one effect \(\beta_j\neq0\)
    regression model better than intercept only

Overall F-test for regression

# mod1 = lm(flipper_length_mm ~ body_mass_g * bill_length_mm, data)
summary(mod1) # same as summary.lm(mod1)
## 
## Call:
## lm(formula = flipper_length_mm ~ body_mass_g * bill_length_mm, 
##     data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.2998  -2.9189   0.0918   2.6606  12.1739 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                 1.983e+02  1.014e+02   1.956   0.0548 .
## body_mass_g                -6.619e-03  2.802e-02  -0.236   0.8140  
## bill_length_mm             -8.085e-01  2.058e+00  -0.393   0.6958  
## body_mass_g:bill_length_mm  3.371e-04  5.647e-04   0.597   0.5527  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.45 on 64 degrees of freedom
## Multiple R-squared:  0.4421, Adjusted R-squared:  0.416 
## F-statistic: 16.91 on 3 and 64 DF,  p-value: 3.389e-08

Overall F-test for regression

summary.lm(mod1)
## 
## Call:
## lm(formula = flipper_length_mm ~ body_mass_g * bill_length_mm, 
##     data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.2998  -2.9189   0.0918   2.6606  12.1739 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                 1.983e+02  1.014e+02   1.956   0.0548 .
## body_mass_g                -6.619e-03  2.802e-02  -0.236   0.8140  
## bill_length_mm             -8.085e-01  2.058e+00  -0.393   0.6958  
## body_mass_g:bill_length_mm  3.371e-04  5.647e-04   0.597   0.5527  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.45 on 64 degrees of freedom
## Multiple R-squared:  0.4421, Adjusted R-squared:  0.416 
## F-statistic: 16.91 on 3 and 64 DF,  p-value: 3.389e-08

Overall F-test for regression

## Multiple R-squared: 0.4421, Adjusted R-squared: 0.416
## F-statistic: 16.91 on 3 and 64 DF, p-value: 3.389e-08

curve(df(x,3,64), from=0, to=17, col=4, lwd=2)
abline(v=qf(p=0.95, df1=3, df2=64), lty=2)
points(16.91, 0, col=2, pch=16, cex=1.5)

Single effect t-test

  • test full regression model   y~x1+x2+x3
    against reduced model       y~ x2+x3
  • H0: full model is only as good as reduced model
           same as \(SSR_f=SSR_r\)
           same as \(\beta_1=0\)
  • F test for \(SSR_f=SSR_r\)    ➔ \(F_{1,n-k-1}\)
    t test for \(\beta_1 / se(\beta_1)=0\)         ➔ \(t_{n-k-1}\)       (identical)
  • If \(P<0.05\) âž” reject H0
  • H1: \(\beta_1\neq0\)
  • ATTN: if x1 involved in interaction x1:x2,
                test for main effect \(\beta_1\) is meaningless

Single effect t-test

## Estimate Std. Error t value Pr(>|t|)
## body_mass_g:bill_length_mm 3.371e-04 5.647e-04 0.597 0.5527

curve(dt(x,64), from=-4, to=4, col=4, lwd=2)
abline(v=qt(p=0.975, df=64), lty=2) # two-sided: alpha/2 = 2.5%
points(0.597, 0, col=2, pch=16, cex=1.5)

Single effect t-test

  • interaction effect in mod1 not significant
  • F-test with anova(mod1, mod2) suggested reduced model mod2
summary(mod2)
## 
## Call:
## lm(formula = flipper_length_mm ~ body_mass_g + bill_length_mm, 
##     data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.2687  -2.9955  -0.0504   2.6024  11.9246 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    1.381e+02  9.859e+00  14.009  < 2e-16 ***
## body_mass_g    1.007e-02  2.009e-03   5.010 4.44e-06 ***
## bill_length_mm 4.122e-01  2.313e-01   1.782   0.0793 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.423 on 65 degrees of freedom
## Multiple R-squared:  0.439,  Adjusted R-squared:  0.4218 
## F-statistic: 25.43 on 2 and 65 DF,  p-value: 6.931e-09

Type I anova() table

multiple regression y~1+x1+...+xk

test model 1 vs model 2
1 y~1 y~1+x1
2 y~1+x1 y~1+x1+x2
3 y~1+x1+x2 y~1+x1+x2+x3
… … …
  • H0: model 1 and model 2 perform equally (\(SSR_1=SSR_2\))
  • If \(P<0.05\) reject H0
  • H1: model 2 performs better than model 1 (\(SSR_1>SSR_2\))

Type I anova() table

multiple regression y~1+x1+...+xk

test model 1 vs model 2
1 y~1 y~1+x1
2 y~1+x1 y~1+x1+x2
3 y~1+x1+x2 y~1+x1+x2+x3
… … …
  • stepwise inclusion of predictors
  • interactions x1:x2 treated as individual predictor
  • order matters for table: y~x1+x2 or y~x2+x1

Type I anova() table

anova(mod1)
## Analysis of Variance Table
## 
## Response: flipper_length_mm
##                            Df  Sum Sq Mean Sq F value    Pr(>F)    
## body_mass_g                 1 1402.68 1402.68 47.2190 3.124e-09 ***
## bill_length_mm              1   93.44   93.44  3.1457   0.08089 .  
## body_mass_g:bill_length_mm  1   10.58   10.58  0.3563   0.55267    
## Residuals                  64 1901.17   29.71                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • 1 + body_mass_g better than 1
  • 1 + body_mass_g + bill_length_mm
    not better than 1 + body_mass_g

Type I anova() table

summary.aov(mod1)
##                            Df Sum Sq Mean Sq F value   Pr(>F)    
## body_mass_g                 1 1402.7  1402.7  47.219 3.12e-09 ***
## bill_length_mm              1   93.4    93.4   3.146   0.0809 .  
## body_mass_g:bill_length_mm  1   10.6    10.6   0.356   0.5527    
## Residuals                  64 1901.2    29.7                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • same output

Type I anova() table

  • Careful with anova(mod1)
  • Type I anova() table usually not recommended in regression with multiple predictors
  • rather use summary()
  • anova(mod1, mod2) F-tests for competing models
  • for ANOVA with categorical predictors
    Type II ANOVA often makes more sense

AIC: Akaike information criterion

\(AIC\)   \(=\)   \(-2\log L\)   \(+\)   \(2k\)           lower=better
                  goodness-of-fit          penalty

  • compare models, lower AIC = better
  • tradeoff between goodness-of-fit and model complexity \(k\)
  • models don’t need to be nested
  • pure value doesn’t say anything (as opposed to \(R^2\)),
    only model difference \(AIC_1-AIC_2\)
  • only compare models on exactly same data,
    watch out for NA in data

BIC: Bayesian information criterion

\(BIC\)   \(=\)   \(-2\log L\)   \(+\)   \(k\log n\)           lower=better
                  goodness-of-fit              penalty

  • \(\log L=\sum_{i=1}^n \log lik(y_i)\) grows with data size \(n\)
  • choose a penalty term that also grows with \(n\)
  • penalizes model complexity \(k\) more strongly than \(AIC\)
  • tends to select for more parsimonious model than \(AIC\)
  • note: BIC does not require Bayesian statistics

AIC and BIC

# mod1 = lm(flipper_length_mm ~ body_mass_g * bill_length_mm, data)
# mod2 = lm(flipper_length_mm ~ body_mass_g , data)
AIC(mod1, mod2)
##      df      AIC
## mod1  5 429.4645
## mod2  4 427.8421

very small difference, still select simpler model

BIC(mod1, mod2)
##      df      BIC
## mod1  5 440.5621
## mod2  4 436.7201

BIC difference larger than AIC difference

Summary regression

  • anova(mod1, mod2): F-test for comparing nested models
  • summary.lm(mod1): F-test for overall model fit and \(R^2\)
  • summary.lm(mod1): t-tests for single parameters
  • summary.aov(mod1): table for stepwise F-tests
    âž” not recommended for regression
  • AIC(mod1, mod2) and BIC(mod1, mod2)

Overview

  1. Hypothesis testing
  2. Regression
  3. ANOVA
  4. ANCOVA
  5. Mixed models

One-way ANOVA

library(palmerpenguins)
data = as.data.frame(penguins)
boxplot(flipper_length_mm ~ species, data=data)

One-way ANOVA

library(ggplot2)
ggplot(data=data, aes(y=flipper_length_mm, x=species)) + 
  geom_boxplot()

One-way ANOVA

response y

one grouping factor A with several levels A1, A2, A3, …

test for difference between group means

H0: all group means identical, \(\bar{y}_1=\bar{y}_2=\bar{y}_3=\ldots\)
âž” identical to overall mean \(\bar{y}\)

F-test / ANOVA for y~1 vs y~A

\(P<0.05\) reject H0

âž” H1: not all group means identical

One-way ANOVA

mod1 = lm(flipper_length_mm ~ species, data=data)
summary.aov(mod1)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## species       2  52473   26237   594.8 <2e-16 ***
## Residuals   339  14953      44                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 2 observations deleted due to missingness

"y~species" is better than "y~1"

One-way ANOVA

# mod1 = lm(flipper_length_mm ~ species, data=data)
summary.lm(mod1)
## 
## Call:
## lm(formula = flipper_length_mm ~ species, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.9536  -4.8235   0.0464   4.8130  20.0464 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      189.9536     0.5405 351.454  < 2e-16 ***
## speciesChinstrap   5.8699     0.9699   6.052 3.79e-09 ***
## speciesGentoo     27.2333     0.8067  33.760  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.642 on 339 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.7782, Adjusted R-squared:  0.7769 
## F-statistic: 594.8 on 2 and 339 DF,  p-value: < 2.2e-16

One-way ANOVA

## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 189.9536 0.5405 351.454 < 2e-16 ***
## speciesChinstrap 5.8699 0.9699 6.052 3.79e-09 ***
## speciesGentoo 27.2333 0.8067 33.760 < 2e-16 ***

  • What are the estimated parameters?
    Group-level means \(\bar{y}_1\), \(\bar{y}_2\), \(\bar{y}_3\), …? ➔ No!
  • lm() uses dummy-coding.
  • (Intercept) group mean of reference level 1 (Adelie)
  • other parameters: differences to reference level
  • t-tests just for differencess to reference level

## F-statistic: 594.8 on 2 and 339 DF, p-value: < 2.2e-16

  • F-test is the important one! Same as in summary.aov()
  • "y~species" is better than "y~1"

One-way ANOVA

But how to get group-level means (predictions)?

Means using model.tables(). Requires aov() object.

aov1 = aov(flipper_length_mm ~ species, data=data)
model.tables(aov1, "means")
## Tables of means
## Grand mean
##          
## 200.9152 
## 
##  species 
##     Adelie Chinstrap Gentoo
##        190     195.8  217.2
## rep    151      68.0  123.0

One-way ANOVA

Contrasts using TukeyHSD(). Requires aov() object.

# aov1 = aov(flipper_length_mm ~ species, data=data)
TukeyHSD(aov1)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = flipper_length_mm ~ species, data = data)
## 
## $species
##                       diff       lwr       upr p adj
## Chinstrap-Adelie  5.869887  3.586583  8.153191     0
## Gentoo-Adelie    27.233349 25.334376 29.132323     0
## Gentoo-Chinstrap 21.363462 19.000841 23.726084     0

TukeyHSD corrects P-values for multiple comparisons.

One-way ANOVA

Plot means and data using sjPlot.

library("sjPlot")
plot_model(mod1, type="pred", show.data = TRUE, jitter=0.5)

One-way ANOVA

emmeans() can predict means and also contrasts,
works with lm() object.

library(emmeans)
emmeans(mod1, ~species)
##  species   emmean    SE  df lower.CL upper.CL
##  Adelie       190 0.540 339      189      191
##  Chinstrap    196 0.805 339      194      197
##  Gentoo       217 0.599 339      216      218
## 
## Confidence level used: 0.95

One-way ANOVA

Contrasts differences using emmeans().

pairs(emmeans(mod1, ~species))
##  contrast           estimate    SE  df t.ratio p.value
##  Adelie - Chinstrap    -5.87 0.970 339  -6.052  <.0001
##  Adelie - Gentoo      -27.23 0.807 339 -33.760  <.0001
##  Chinstrap - Gentoo   -21.36 1.004 339 -21.286  <.0001
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

One-way ANOVA

Plot means using emmeans().

plot(emmeans(mod1, ~species), horizontal=FALSE)

One-way ANOVA

Plot means using emmip().

emmip(mod1, ~species, CIs=TRUE) +
  geom_jitter( aes(x=species, y=flipper_length_mm), 
               data=data, width=0.1, alpha=0.2, color="red")

One-way ANOVA - conclusions

  • ANOVA is just a linear model
  • F-test for difference in group-means
  • y~A vs y~1
  • lm() and aov() use dummy-coding.
    estimated parameters are not group-means
  • predicted group-means (emmeans) are model-dependent!
    important in two-way ANOVA and GLM
  • post-hoc tests / contrast:
    test for differences in certain groups (t-tests)
    P-values are corrected for multiple comparisons (e.g. Tukey)
  • with a GLM, use anova( , test="Chi") for a Chi-squared likelihood-ratio-test (LRT) instead

Overview

  1. Hypothesis testing
  2. Regression
  3. ANOVA
  4. ANCOVA
  5. Mixed models

Two-way ANOVA

data = subset(data, !is.na(sex))
ggplot(data = data, aes(y=flipper_length_mm, x=species, colour=sex)) + 
  geom_boxplot()

Two-way ANOVA   y~A+B

mod1 = lm(flipper_length_mm ~ species + sex, data=data)
summary.lm(mod1)
## 
## Call:
## lm(formula = flipper_length_mm ~ species + sex, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5278  -3.7239   0.2761   3.4722  16.4722 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      186.6777     0.5684 328.437  < 2e-16 ***
## speciesChinstrap   5.7208     0.8407   6.805 4.79e-11 ***
## speciesGentoo     27.0462     0.7072  38.243  < 2e-16 ***
## sexmale            6.8502     0.6276  10.914  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.726 on 329 degrees of freedom
## Multiple R-squared:  0.8346, Adjusted R-squared:  0.8331 
## F-statistic: 553.4 on 3 and 329 DF,  p-value: < 2.2e-16

Two-way ANOVA   y~A+B

mod1 = lm(flipper_length_mm ~ species + sex, data=data)

estimated group-means:

B1 B2
A1 (intercept) (intercept) + b2
A2 (intercept) + a2 (intercept) + a2 + b2 + a2:b2
  • effect of sex independent of species
  • effect of species independent of sex
  • Df < level combinations

Two-way ANOVA   y~A*B

mod2 = lm(flipper_length_mm ~ species * sex, data=data)
summary.lm(mod2)
## 
## Call:
## lm(formula = flipper_length_mm ~ species * sex, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.7945  -3.4110   0.0882   3.4590  17.5890 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              187.7945     0.6619 283.721  < 2e-16 ***
## speciesChinstrap           3.9408     1.1742   3.356 0.000884 ***
## speciesGentoo             24.9124     0.9947  25.044  < 2e-16 ***
## sexmale                    4.6164     0.9361   4.932  1.3e-06 ***
## speciesChinstrap:sexmale   3.5600     1.6606   2.144 0.032782 *  
## speciesGentoo:sexmale      4.2176     1.3971   3.019 0.002737 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.655 on 327 degrees of freedom
## Multiple R-squared:  0.8396, Adjusted R-squared:  0.8372 
## F-statistic: 342.4 on 5 and 327 DF,  p-value: < 2.2e-16

Two-way ANOVA   y~A*B

mod2 = lm(flipper_length_mm ~ species * sex, data=data)

estimated group-means:

B1 B2
A1 (intercept) (intercept) + b2
A2 (intercept) + a2 (intercept) + a2 + b2 + a2:b2
  • effect of sex changes with species
  • effect of species changes with sex
  • Df = level combinations
  • fits independent means for all level combinations

Type 1 anova() table

two-way ANOVA w. interaction y~A*B, or y~1+A+B+A:B

test model 1 vs model 2
A y~1 y~1+A
B y~1+A y~1+A+B
A:B y~1+A+B y~1+A+B+A:B
  • includes predictors stepwise
  • tests against previous model
  • order matters for the table

Type 1 anova() table

# mod2 = lm(flipper_length_mm ~ species * sex, data=data)
summary.aov(mod2) # same as anova(mod2)
##              Df Sum Sq Mean Sq F value  Pr(>F)    
## species       2  50526   25263 789.912 < 2e-16 ***
## sex           1   3906    3906 122.119 < 2e-16 ***
## species:sex   2    329     165   5.144 0.00631 ** 
## Residuals   327  10458      32                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Type 2 Anova() table

two-way ANOVA w. interaction y~A*B, or y~1+A+B+A:B

test reduced model vs full model
A y~1+ B y~1+A+B+A:B
B y~1+A y~1+A+B+A:B
A:B y~1+A+B y~1+A+B+A:B
  • independently removes predictors completely
    (with its higher-order interactions)
  • tests against full model
  • order does not matter for the table

Type 2 Anova() table

# mod2 = lm(flipper_length_mm ~ species * sex, data=data)
library(car)
Anova(mod2) # same as default: Anova(mod2, type=2)
## Anova Table (Type II tests)
## 
## Response: flipper_length_mm
##             Sum Sq  Df  F value    Pr(>F)    
## species      50185   2 784.5829 < 2.2e-16 ***
## sex           3906   1 122.1189 < 2.2e-16 ***
## species:sex    329   2   5.1442  0.006314 ** 
## Residuals    10458 327                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Two-way ANOVA

predicted means

# mod2 = lm(flipper_length_mm ~ species * sex, data=data)
library(emmeans)
emmeans(mod2, ~species*sex)
##  species   sex    emmean    SE  df lower.CL upper.CL
##  Adelie    female    188 0.662 327      186      189
##  Chinstrap female    192 0.970 327      190      194
##  Gentoo    female    213 0.743 327      211      214
##  Adelie    male      192 0.662 327      191      194
##  Chinstrap male      200 0.970 327      198      202
##  Gentoo    male      222 0.724 327      220      223
## 
## Confidence level used: 0.95
# emmeans(mod2, ~species|sex) # same, just ordered
# emmeans(mod2, ~sex|species) # same, just ordered the other way

Two-way ANOVA

marginal means

# mod2 = lm(flipper_length_mm ~ species * sex, data=data)
emmeans(mod2, ~species)
## NOTE: Results may be misleading due to involvement in interactions
##  species   emmean     SE  df lower.CL upper.CL
##  Adelie     190.1 0.4680 327    189.2    191.0
##  Chinstrap  195.8 0.6858 327    194.5    197.2
##  Gentoo     217.1 0.5186 327    216.1    218.1
## 
## Results are averaged over the levels of: sex 
## Confidence level used: 0.95
emmeans(mod2, ~sex)
## NOTE: Results may be misleading due to involvement in interactions
##  sex    emmean     SE  df lower.CL upper.CL
##  female  197.4 0.4631 327    196.5    198.3
##  male    204.6 0.4598 327    203.7    205.5
## 
## Results are averaged over the levels of: species 
## Confidence level used: 0.95

Two-way ANOVA

all contrasts, as list

pairs( emmeans(mod2, ~species*sex) )
##  contrast                          estimate    SE  df t.ratio p.value
##  Adelie female - Chinstrap female    -3.941 1.174 327  -3.356  0.0113
##  Adelie female - Gentoo female      -24.912 0.995 327 -25.044  <.0001
##  Adelie female - Adelie male         -4.616 0.936 327  -4.932  <.0001
##  Adelie female - Chinstrap male     -12.117 1.174 327 -10.320  <.0001
##  Adelie female - Gentoo male        -33.746 0.981 327 -34.399  <.0001
##  Chinstrap female - Gentoo female   -20.972 1.221 327 -17.169  <.0001
##  Chinstrap female - Adelie male      -0.676 1.174 327  -0.575  0.9926
##  Chinstrap female - Chinstrap male   -8.176 1.372 327  -5.961  <.0001
##  Chinstrap female - Gentoo male     -29.806 1.210 327 -24.626  <.0001
##  Gentoo female - Adelie male         20.296 0.995 327  20.403  <.0001
##  Gentoo female - Chinstrap male      12.795 1.221 327  10.475  <.0001
##  Gentoo female - Gentoo male         -8.834 1.037 327  -8.518  <.0001
##  Adelie male - Chinstrap male        -7.501 1.174 327  -6.388  <.0001
##  Adelie male - Gentoo male          -29.130 0.981 327 -29.694  <.0001
##  Chinstrap male - Gentoo male       -21.629 1.210 327 -17.870  <.0001
## 
## P value adjustment: tukey method for comparing a family of 6 estimates

Two-way ANOVA

all contrasts, as table

pwpm( emmeans(mod2, ~species*sex) )
##                  Adelie female Chinstrap female Gentoo female Adelie male
## Adelie female            [188]           0.0113        <.0001      <.0001
## Chinstrap female        -3.941            [192]        <.0001      0.9926
## Gentoo female          -24.912          -20.972         [213]      <.0001
## Adelie male             -4.616           -0.676        20.296       [192]
## Chinstrap male         -12.117           -8.176        12.795      -7.501
## Gentoo male            -33.746          -29.806        -8.834     -29.130
##                  Chinstrap male Gentoo male
## Adelie female            <.0001      <.0001
## Chinstrap female         <.0001      <.0001
## Gentoo female            <.0001      <.0001
## Adelie male              <.0001      <.0001
## Chinstrap male            [200]      <.0001
## Gentoo male             -21.629       [222]
## 
## Row and column labels: species:sex
## Upper triangle: P values   adjust = "tukey"
## Diagonal: [Estimates] (emmean) 
## Lower triangle: Comparisons (estimate)   earlier vs. later

Two-way ANOVA

in-group contrasts only

pairs( emmeans(mod2, ~species|sex) )
## sex = female:
##  contrast           estimate    SE  df t.ratio p.value
##  Adelie - Chinstrap    -3.94 1.174 327  -3.356  0.0025
##  Adelie - Gentoo      -24.91 0.995 327 -25.044  <.0001
##  Chinstrap - Gentoo   -20.97 1.221 327 -17.169  <.0001
## 
## sex = male:
##  contrast           estimate    SE  df t.ratio p.value
##  Adelie - Chinstrap    -7.50 1.174 327  -6.388  <.0001
##  Adelie - Gentoo      -29.13 0.981 327 -29.694  <.0001
##  Chinstrap - Gentoo   -21.63 1.210 327 -17.870  <.0001
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

Two-way ANOVA

in-group contrasts only

pairs( emmeans(mod2, ~sex|species) )
## species = Adelie:
##  contrast      estimate    SE  df t.ratio p.value
##  female - male    -4.62 0.936 327  -4.932  <.0001
## 
## species = Chinstrap:
##  contrast      estimate    SE  df t.ratio p.value
##  female - male    -8.18 1.372 327  -5.961  <.0001
## 
## species = Gentoo:
##  contrast      estimate    SE  df t.ratio p.value
##  female - male    -8.83 1.037 327  -8.518  <.0001

Two-way ANOVA

marginal contrasts

pairs( emmeans(mod2, ~sex) )
##  contrast      estimate    SE  df t.ratio p.value
##  female - male    -7.21 0.653 327 -11.046  <.0001
## 
## Results are averaged over the levels of: species
pairs( emmeans(mod2, ~species) )
##  contrast           estimate    SE  df t.ratio p.value
##  Adelie - Chinstrap    -5.72 0.830 327  -6.890  <.0001
##  Adelie - Gentoo      -27.02 0.699 327 -38.681  <.0001
##  Chinstrap - Gentoo   -21.30 0.860 327 -24.774  <.0001
## 
## Results are averaged over the levels of: sex 
## P value adjustment: tukey method for comparing a family of 3 estimates

Two-way ANOVA

marginal contrasts

mod2 |> emmeans(~sex) |> pairs()
##  contrast      estimate    SE  df t.ratio p.value
##  female - male    -7.21 0.653 327 -11.046  <.0001
## 
## Results are averaged over the levels of: species
mod2 |> emmeans(~species) |> pairs()
##  contrast           estimate    SE  df t.ratio p.value
##  Adelie - Chinstrap    -5.72 0.830 327  -6.890  <.0001
##  Adelie - Gentoo      -27.02 0.699 327 -38.681  <.0001
##  Chinstrap - Gentoo   -21.30 0.860 327 -24.774  <.0001
## 
## Results are averaged over the levels of: sex 
## P value adjustment: tukey method for comparing a family of 3 estimates

Two-way ANOVA

emmip(mod2, ~species|sex, CIs=TRUE) +
  geom_jitter( aes(x=species, y=flipper_length_mm), 
             data=data, width=0.1, alpha=0.2, color="red")

Two-way ANOVA

emmip(mod2, ~sex|species, CIs=TRUE) +
  geom_jitter( aes(x=sex, y=flipper_length_mm), 
             data=data, width=0.1, alpha=0.2, color="red")

Two-way ANOVA - conclusions

  • test if means change with factor levels
  • y~A+B vs y~A*B
  • Type 2 Anova() table for multiple predictors or interactions
  • model-based group means with emmeans
    e.g. y~A+B and y~A*B have different means
  • post-hoc tests / contrast (e.g. Tukey)
  • with a GLM, Anova() automatically chooses a
    Chi-squared likelihood-ratio-test (LRT) instead

Overview

  1. Hypothesis testing
  2. Regression
  3. ANOVA
  4. ANCOVA
  5. Mixed models

ANCOVA

ggplot( aes(x=body_mass_g, y=flipper_length_mm, color=species), data=data) + 
  geom_point()

ANCOVA   y~A+x

data$mass_z = as.numeric( scale(data$body_mass_g) )
mod1 = lm(flipper_length_mm ~ mass_z + species, data=data)
summary.lm(mod1)
## 
## Call:
## lm(formula = flipper_length_mm ~ mass_z + species, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.5418  -3.1804   0.0983   3.3295  17.3954 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      194.3677     0.5520 352.092  < 2e-16 ***
## mass_z             6.8561     0.5200  13.186  < 2e-16 ***
## speciesChinstrap   5.4915     0.7938   6.918  2.4e-11 ***
## speciesGentoo     15.3289     1.1167  13.727  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.405 on 329 degrees of freedom
## Multiple R-squared:  0.8526, Adjusted R-squared:  0.8513 
## F-statistic: 634.4 on 3 and 329 DF,  p-value: < 2.2e-16

ANCOVA   y~A+x

group intercept slope
A1 (intercept) x
A2 (intercept) + a2 x + x:a2
A3 (intercept) + a3 x + x:a3
  • individual intercepts per group (dummy-coded)
  • slope identical across groups
  • test for difference in group-means while correcting for covariate effect: y~A+x vs y~x

ANCOVA   y~A+x

library("sjPlot")
plot_model(mod1, 
           type="pred", 
           terms=c("mass_z", "species"), 
           show.data = TRUE) + ylim(170, 240)

ANCOVA   y~A*x

# data$mass_z = as.numeric( scale(data$body_mass_g) )
mod2 = lm(flipper_length_mm ~ mass_z * species, data=data)
summary.lm(mod2)
## 
## Call:
## lm(formula = flipper_length_mm ~ mass_z * species, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.4296  -3.2343   0.1668   3.2517  17.9549 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             193.4139     0.6572 294.302  < 2e-16 ***
## mass_z                    5.3229     0.7804   6.820 4.40e-11 ***
## speciesChinstrap          8.0523     1.2262   6.567 2.02e-10 ***
## speciesGentoo            15.5511     1.1956  13.007  < 2e-16 ***
## mass_z:speciesChinstrap   4.2633     1.5767   2.704  0.00721 ** 
## mass_z:speciesGentoo      2.1986     1.1113   1.978  0.04872 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.352 on 327 degrees of freedom
## Multiple R-squared:  0.8564, Adjusted R-squared:  0.8542 
## F-statistic: 389.9 on 5 and 327 DF,  p-value: < 2.2e-16

ANCOVA   y~A*x

group intercept slope
A1 (intercept) x
A2 (intercept) + a2 x + x:a2
A3 (intercept) + a3 x + x:a3
  • individual intercepts per group (dummy-coded)
  • individual slopes per group (dummy-coded)
  • test for difference in slopes: y~A*x vs y~A+x

ANCOVA   y~A*x

# library("sjPlot")
plot_model(mod2, 
           type="pred", 
           terms=c("mass_z", "species"), 
           show.data = TRUE) + ylim(170, 240)

Type 2 Anova() table

# mod2 = lm(flipper_length_mm ~ mass_z * species, data=data)
Anova(mod2) 
## Anova Table (Type II tests)
## 
## Response: flipper_length_mm
##                Sum Sq  Df F value  Pr(>F)    
## mass_z         5080.0   1 177.319 < 2e-16 ***
## species        5903.2   2 103.028 < 2e-16 ***
## mass_z:species  244.6   2   4.269 0.01478 *  
## Residuals      9368.2 327                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANCOVA

estimated group-level slopes

# mod2 = lm(flipper_length_mm ~ mass_z * species, data=data)
library("emmeans")
emtrends(mod2, ~species, var="mass_z")
##  species   mass_z.trend    SE  df lower.CL upper.CL
##  Adelie            5.32 0.780 327     3.79     6.86
##  Chinstrap         9.59 1.370 327     6.89    12.28
##  Gentoo            7.52 0.791 327     5.97     9.08
## 
## Confidence level used: 0.95

ANCOVA

estimated contrasts for group-level slopes

# mod2 = lm(flipper_length_mm ~ mass_z * species, data=data)
pairs( emtrends(mod2, ~species, var="mass_z") )
##  contrast           estimate   SE  df t.ratio p.value
##  Adelie - Chinstrap    -4.26 1.58 327  -2.704  0.0197
##  Adelie - Gentoo       -2.20 1.11 327  -1.978  0.1192
##  Chinstrap - Gentoo     2.06 1.58 327   1.305  0.3933
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

Overview

  1. Hypothesis testing
  2. Regression
  3. ANOVA
  4. ANCOVA
  5. Mixed models

Mixed effects models

again, continuous predictor(s) x, categorical predictor A

  • fixed effects only model:

     y~x+A separate intercepts model (ANCOVA)
     y~x*A separate intercepts & slopes model (ANCOVA)

  • mixed effects model:

     y~x+(1|A) random intercepts model (LMM)
     y~x+(1+x|A) random intercepts & slopes model (LMM)

âž” How to find best model (fixed and random effects structure) with AIC / BIC / F-tests etc?

LMM protocol

library(lme4)
lmer( y ~ x + (1+x|A), data=data, REML=TRUE )

Zuur et al.: Mixed effects models (2009)

  1. Start with maximal fixed effects structure
  2. Find optimal random effects structure using REML=TRUE
  3. Find optimal fixed effects structure using REML=FALSE
  4. Fit best model again using REML=TRUE

LMM protocol - step 1

ANCOVA model from last section, add sex as random variable

lmmax1 = lmer( flipper_length_mm ~ mass_z * species + 
                 (1+mass_z|sex), data=data )

Fixed effects have maximal structure.

LMM protocol - step 2

ranova() tests if some random effects can be dropped.
Uses Chi-sq-test for likelihood-ratio instead of F-test for residuals

library("lmerTest")
ranova(lmmax1)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## flipper_length_mm ~ mass_z + species + (1 + mass_z | sex) + mass_z:species
##                              npar  logLik    AIC    LRT Df Pr(>Chisq)
## <none>                         10 -1021.1 2062.1                     
## mass_z in (1 + mass_z | sex)    8 -1021.8 2059.6 1.4519  2     0.4839

H0: Full model is equal to reduced model.

\(P>0.05\) âž” drop (mass_z|sex) term (random slopes)

LMM protocol - step 2

lmmax2 = lmer( flipper_length_mm ~ mass_z * species + 
                 (1|sex), data=data )
ranova(lmmax2)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## flipper_length_mm ~ mass_z + species + (1 | sex) + mass_z:species
##           npar  logLik    AIC    LRT Df Pr(>Chisq)  
## <none>       8 -1021.8 2059.6                       
## (1 | sex)    7 -1024.8 2063.5 5.9104  1    0.01505 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

H0: Full model is equal to reduced model.

\(P<0.05\) âž” keep (1|sex) term (random intercepts)

LMM protocol - step 3

anova() tests if some fixed effects can be dropped.
Uses Chi-sq-test for likelihood-ratio instead of F-test for residuals

Fit models with REML=FALSE !

lmm1 = lmer( flipper_length_mm ~ mass_z * species + (1|sex), data=data, REML=FALSE )
lmm2 = lmer( flipper_length_mm ~ mass_z + species + (1|sex), data=data, REML=FALSE )
anova(lmm1, lmm2)
## Data: data
## Models:
## lmm2: flipper_length_mm ~ mass_z + species + (1 | sex)
## lmm1: flipper_length_mm ~ mass_z * species + (1 | sex)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## lmm2    6 2074.5 2097.3 -1031.2   2062.5                        
## lmm1    8 2069.3 2099.7 -1026.6   2053.3 9.2154  2   0.009975 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

H0: Models perform equally well.

\(P<0.05\) âž” Keep full model.

LMM protocol - step 4

refit model with REML=TRUE (default in lmer())

lmm1 = lmer( flipper_length_mm ~ mass_z * species + (1|sex), data=data)
coef(summary(lmm1)) |> round(3)
##                         Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)              192.487      1.413   1.442 136.196    0.001
## mass_z                     3.833      0.929 233.914   4.125    0.000
## speciesChinstrap           8.211      1.212 326.158   6.773    0.000
## speciesGentoo             18.085      1.474 212.634  12.265    0.000
## mass_z:speciesChinstrap    4.448      1.559 326.131   2.854    0.005
## mass_z:speciesGentoo       2.200      1.098 326.000   2.004    0.046

LMM protocol - step 4

fixef(lmm1) 
##             (Intercept)                  mass_z        speciesChinstrap 
##              192.487198                3.833164                8.210829 
##           speciesGentoo mass_z:speciesChinstrap    mass_z:speciesGentoo 
##               18.085338                4.448020                2.200067
ranef(lmm1) 
## $sex
##        (Intercept)
## female   -1.145652
## male      1.145652
## 
## with conditional variances for "sex"

LMM protocol - step 4

library("sjPlot")
plot_model(lmm1, 
           type="pred", 
           terms=c("mass_z", "species"), 
           show.data = TRUE) 

LMM - automatize steps 2 & 3

lmmax1 = lmer( flipper_length_mm ~ mass_z * species + (1+mass_z|sex), data=data )
step(lmmax1)
## Backward reduced random-effect table:
## 
##                              Eliminated npar  logLik    AIC    LRT Df
## <none>                                    10 -1021.1 2062.1          
## mass_z in (1 + mass_z | sex)          1    8 -1021.8 2059.6 1.4519  2
## (1 | sex)                             0    7 -1024.8 2063.5 5.9104  1
##                              Pr(>Chisq)  
## <none>                                   
## mass_z in (1 + mass_z | sex)    0.48386  
## (1 | sex)                       0.01505 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                Eliminated Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## mass_z:species          0 260.24  130.12     2 326.07  4.6564 0.01014 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## flipper_length_mm ~ mass_z + species + (1 | sex) + mass_z:species

LMM - emmeans

library("emmeans")
emtrends(lmm1, ~species, var="mass_z")
##  species   mass_z.trend    SE  df lower.CL upper.CL
##  Adelie            3.83 0.997 234     1.87     5.80
##  Chinstrap         8.28 1.462 318     5.40    11.16
##  Gentoo            6.03 1.005 237     4.05     8.01
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
pairs( emtrends(lmm1, ~species, var="mass_z") )
##  contrast           estimate   SE  df t.ratio p.value
##  Adelie - Chinstrap    -4.45 1.56 326  -2.853  0.0128
##  Adelie - Gentoo       -2.20 1.10 326  -2.004  0.1127
##  Chinstrap - Gentoo     2.25 1.56 326   1.437  0.3232
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

Consultation hour

  • Theory in Biodiversity Science group (TiBS) provides consultation in theory, maths and statistics.
  • each Wednesday after the seminar
  • write me an email with a brief description of your problem
    benjamin.rosenbaum@idiv.de

That’s it!

Artwork by @allison_horst